#ifndef __jacobi_fe_hh__
#define __jacobi_fe_hh__

#include "modal_fe.hh"
#include "utils.hh"
#include "femlib_exceptions.hh"

#include "fe_evaluators.hh"

// ALL Jacobi polynomial type 1-D basis functions: nth Order elements
template <unsigned int ORDER, typename T=double>
class JacobiSegmentFiniteElement : public ModalFiniteElement<1, fastfem::Segment, ORDER, fastfem::Jacobi, T> 
{
public:
   JacobiSegmentFiniteElement(double p_alpha=0.0, double p_beta=0.0) : ModalFiniteElement<1, fastfem::Segment, ORDER, fastfem::Jacobi, T> ()
   { }

   virtual ~JacobiSegmentFiniteElement()
   { }

  /** virtual function which evaluates the values of all
       shape functions at a given point ip */
  virtual T evaluate_basis ( unsigned int , const Point<SPATIAL_DIMENSION, T>& ) ;

  /** virtual function which evaluates the values of all
       shape functions at all point ips */
  virtual void evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<T>& basis_vals ) ;

  /** virtual function which evaluates the values of all
       shape functions at a given point ip */
  template <typename U>
  static U evaluate_basis ( unsigned int , const Point<SPATIAL_DIMENSION, U>& ) ;

  /** virtual function which evaluates the values of all
       shape functions at all point ips */
  template <typename U>
  static void evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, U>& ips, std::vector<U>& basis_vals ) ;

  virtual void evaluate_gradient ( unsigned int , const Point<SPATIAL_DIMENSION, T>&, Point<SPATIAL_DIMENSION, T>& pt) ;

  virtual void evaluate_gradient_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<Point<SPATIAL_DIMENSION, T> >& grad_vals ) ;

  virtual void evaluate_hessian ( unsigned int , const Point<SPATIAL_DIMENSION, T>&, PointList<SPATIAL_DIMENSION, T>& pt) ;

  virtual void evaluate_hessian_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<PointList<SPATIAL_DIMENSION, T> >& grad_vals ) ;

  void reference_element_coords(std::vector<Point<SPATIAL_DIMENSION, T> >& nodes) ;  

  virtual int test() ;

  static const std::vector<Point<SPATIAL_DIMENSION, T> >& get_vertices() ;
  
protected:
  const double alpha, beta;   // Jacobi-polynomial parameters
  // Legendre: alpha = beta = 0.0

  static std::vector<Point<SPATIAL_DIMENSION, T> > vertices;
  static std::vector<T> weights;
  
  static void compute_vertices(std::vector<Point<SPATIAL_DIMENSION, T> > &) ;
};


template <unsigned int ORDER, typename T>
struct FEMTraits<JacobiSegmentFiniteElement<ORDER, T> >
{
  
  enum { FE_Dimension = 1 } ;
  enum { FE_Family = fastfem::Jacobi } ;
  enum { FE_Order = ORDER } ;
  enum { N_Modes = ORDER+1 } ;

  typedef T BaseType;
  typedef JacobiSegmentFiniteElement<ORDER, T> FEType;
  typedef NodalFiniteElement<1, fastfem::Segment, ORDER, fastfem::Jacobi> ParentFEType;

};

// All inline functions

/** virtual function which evaluates the values of all
      shape functions at a given point ip */
template <unsigned int ORDER, typename T>
inline
T JacobiSegmentFiniteElement<ORDER, T>::evaluate_basis (  unsigned int ib, const Point<SPATIAL_DIMENSION, T>& px)
{
  femlib_assert(ib < ORDER+1, "Invalid basis index request");

  T basis_val = 0.0;

  Polylib::jacobfd(numPoints, z.data(), mode, NULL, ib, alpha, beta);
  // normalise
  basis_val *= std::sqrt(0.5*(2.0*ib+1.0));

  return basis_val ;
}

/** virtual function which evaluates the values of all
      shape functions at all point ips */
template<unsigned int ORDER, typename T>
inline
void JacobiSegmentFiniteElement<ORDER, T>::evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, T>& ips, std::vector<T>& basis_vals )
{
  basis_vals.resize((ORDER+1)*ips.n_points, 1.0);

  for (unsigned int k = 0, offset=0; k < ips.n_points; k++, offset+=(ORDER+1) )
    for (unsigned int ib = 0; ib < ORDER+1; ib++ )
      basis_vals[ib+offset] = femlib::pow<T>(ips.x(k), ib) ;
  
  return ;
}


template <unsigned int ORDER, typename T>
inline
void JacobiSegmentFiniteElement<ORDER, T>::evaluate_gradient ( unsigned int ib, 
                                                                 const Point<SPATIAL_DIMENSION, T>& px, 
                                                                 Point<SPATIAL_DIMENSION, T>& grad)
{
  femlib_assert(ib < ORDER+1, "Invalid basis index request");
    
  if (ib == 0)
    grad[0] = 0.0 ;
  else if (ib == 1)
    grad[0] = 1.0 ;
  else
    grad[0] = ib * femlib::pow<T>(px.x(), ib-1) ;

  return ;
}

template <unsigned int ORDER, typename T>
inline
void JacobiSegmentFiniteElement<ORDER, T>::evaluate_gradient_all ( const PointList<SPATIAL_DIMENSION, T>& ips, 
                                                                      std::vector<Point<SPATIAL_DIMENSION, T> >& grad_vals )
{
  grad_vals.resize((ORDER+1)*ips.n_points);

  for (unsigned int l = 0, offset=0; l < ips.n_points; l++, offset+=(ORDER+1) )
  {
    for (unsigned int ib = 0; ib < ORDER+1; ib++ )
    {
      if (ib == 0)
        grad_vals[offset+ib][0] = 0.0 ;
      else if (ib == 1)
        grad_vals[offset+ib][0] = 1.0 ;
      else
        grad_vals[offset+ib][0] = ib * femlib::pow<T>(ips.x(l), ib-1) ;
    }
  }
  
  return ;
}


template <unsigned int ORDER, typename T>
inline
void JacobiSegmentFiniteElement<ORDER, T>::evaluate_hessian ( unsigned int ib, 
                                                              const Point<SPATIAL_DIMENSION, T>& px, 
                                                              PointList<SPATIAL_DIMENSION, T>& hess)
{
  femlib_assert(ib < ORDER+1, "Invalid basis index request");
  
  if (ib < 2)
    hess.loc[0] = 0.0 ;
  else if (ib == 2)
    hess.loc[0] = ib*(ib-1) ;
  else
    hess.loc[0] = ib*(ib-1)*femlib::pow<T>(px.x(), ib-2) ;

  return ;
}


template <unsigned int ORDER, typename T>
inline
void JacobiSegmentFiniteElement<ORDER, T>::evaluate_hessian_all ( const PointList<SPATIAL_DIMENSION, T>& ips, 
                                                                  std::vector<PointList<SPATIAL_DIMENSION, T> >& hess_vals )
{
  hess_vals.resize((ORDER+1)*ips.n_points);
  for (unsigned int l = 0; l < ips.n_points; l++ )
  {
    for (unsigned int ib = 0; ib < ORDER+1; ib++ )
    {
      T* pval = hess_vals[l*(ORDER+1)+ib].get(0);
        
      if (ib < 2)
        pval[0] = 0.0 ;
      else if (ib == 2)
        pval[0] = ib*(ib-1) ;
      else
        pval[0] = ib*(ib-1)*femlib::pow<T>(ips.x(l), ib-2) ;
    }
  }

  return ;
}
  


/** virtual function which evaluates the values of all
      shape functions at a given point ip */
template <unsigned int ORDER, typename T>
template <typename U>
U JacobiSegmentFiniteElement<ORDER, T>::evaluate_basis (  unsigned int ib, const Point<SPATIAL_DIMENSION, U>& px)
{
  femlib_assert(ib < ORDER+1, "Invalid basis index request");

  U exponent = static_cast<T>(ib);
  return femlib::pow<U>(px.x(), exponent) ;
}

/** virtual function which evaluates the values of all
      shape functions at all point ips */
template<unsigned int ORDER, typename T>
template <typename U>
void JacobiSegmentFiniteElement<ORDER, T>::evaluate_basis_all ( const PointList<SPATIAL_DIMENSION, U>& ips, std::vector<U>& basis_vals )
{
  basis_vals.resize((ORDER+1)*ips.n_points, 1.0);
  
  for (unsigned int k = 0, offset=0; k < ips.n_points; k++, offset+=(ORDER+1) )
    for (unsigned int ib = 0; ib < ORDER+1; ib++ )
    {
      U exponent = static_cast<T>(ib);

      basis_vals[ib+offset] = femlib::pow<U>(ips.x(k), exponent) ;
    }
  
  return ;
}


template <unsigned int ORDER, typename T>
inline
void JacobiSegmentFiniteElement<ORDER, T>::reference_element_coords(std::vector<Point<SPATIAL_DIMENSION, T> >& nodes)
{
  nodes.resize(2);

  nodes[0][0] = -1.0;
  nodes[1][0] = 1.0 ;

  return ;
}


template <unsigned int ORDER, typename T>
int JacobiSegmentFiniteElement<ORDER, T>::test ( )
{
  int i, j;
  int test_num = 5;
  int nmodes = ORDER+1;

  std::vector<Point<SPATIAL_DIMENSION, T> > nodes ;
  this->reference_element_coords(nodes);

  cout << "\n\nBASIS_MONOMIAL_SEGMENT_TEST: Order=" << ORDER << "\n";
  cout << "  Verify basis functions for 1-D elements.\n";

  PointList<SPATIAL_DIMENSION, T> ips(test_num);

  for ( j = 0; j < test_num; j++ )
    ips.x(j) = generate_random();

  cout << "\n   Number of random evaluation nodes:= " << test_num ;
  for (i=0; i < test_num; ++i)
    cout << "\n        nodes[" << i << "] = [" << ips.x(i)  << ", " << ips.y(i)  << ", " << ips.z(i) << "]" ;
  cout << endl << endl ;

  
  cout << "\n  The basis function values at random nodes (Vandermonde matrix)\n\n";

  for ( j = 0; j < test_num; j++ )
  {
    Point<SPATIAL_DIMENSION, T> pt = ips.get_point(j);

    cout << "  ";
    for ( i = 0; i < nmodes; i++ )
    {
      T phi = this->evaluate_basis ( i, pt );
      cout << setw(12) << phi;
    }
    cout << "\n";
  }
  
  cout << "\n  The X derivatives should sum to 0.\n\n";
  cout << "  \t\tdPhi_i/dX_j ... \t\tsum(dPhi/dX):\n\n";

  // evaluate_gradient (  unsigned int ib, const Point<SPATIAL_DIMENSION>& px, Point<SPATIAL_DIMENSION>& grad)
  Point<SPATIAL_DIMENSION, T> grad ;
  for ( j = 0; j < test_num; j++ )
  {
    Point<SPATIAL_DIMENSION, T> pt = ips.get_point(j);

    T sum_x = 0.0;
    cout << "  ";
    for ( i = 0; i < nmodes; i++ )
    {
      this->evaluate_gradient(i, pt, grad);
      cout << setw(12) << grad[0];
      sum_x += grad[0];
    }
    cout << setw(15) << sum_x << "\n";
  }
  
  cout << "\n  The X double-derivatives should sum to 0.\n\n";
  cout << "  \t\td^2Phi_i/dX^2_j ... \t\tsum(d^2Phi/dX^2):\n\n";

  // evaluate_gradient (  unsigned int ib, const Point<SPATIAL_DIMENSION>& px, Point<SPATIAL_DIMENSION>& grad)
  PointList<SPATIAL_DIMENSION, T> hess(SPATIAL_DIMENSION) ;
  for ( j = 0; j < test_num; j++ )
  {
    Point<SPATIAL_DIMENSION, T> pt = ips.get_point(j);

    T sum_x = 0.0;
    cout << "  ";
    for ( i = 0; i < nmodes; i++ )
    {
      this->evaluate_hessian(i, pt, hess);
      cout << setw(12) << hess.x(0);
      sum_x += hess.x(0);
    }
    cout << setw(15) << sum_x << "\n";
  }
  

  cout << "\n==   The basis function values at ANY point P should sum to 1   ==\n";
  cout << "   ---------------------P---------------------------     PHI_SUM\n";

  // evaluate_basis_all ( const PointList<SPATIAL_DIMENSION>& ips, std::vector<double>& basis_vals )
  std::vector<T> basis_vals;
  this->evaluate_basis_all/*_unrolled*/ ( ips, basis_vals );

  for ( j = 0; j < test_num; j++ )
  {
    T phi_sum = 0.0;
    for (i=0; i < nmodes; ++i)
    {
      cout << "  " << setw(12) << basis_vals[i+j*nmodes] ;
      phi_sum += basis_vals[i+j*nmodes] ;
    }
    cout << "  " << setw(12) << phi_sum << "\n";
  }

  basis_vals.clear();
  
  cout << "\n==   The basis derivative values at ANY point P should sum to 0   ==\n";
  cout << "   ---------------------P-----------------------     PHI_DERIVATIVE_SUM\n";

  std::vector<Point<SPATIAL_DIMENSION, T> > grad_vals;
  this->evaluate_gradient_all(ips, grad_vals);

  for ( j = 0; j < test_num; j++ )
  {
    T dphi_sum = 0.0 ;
    for (i=0; i < nmodes; ++i)
    {
      cout << "  " << setw(12) << grad_vals[j*nmodes+i][0] ;
      dphi_sum += grad_vals[j*nmodes+i][0] ;
    }
    cout << "  " << setw(12) << dphi_sum << "\n";
  }

  grad_vals.clear();

  cout << "\n==   The basis second-derivative values at ANY point P should sum to 0   ==\n";
  cout << "   ---------------------P-----------------------     DPHI_DERIVATIVE_SUM\n";

  std::vector<PointList<SPATIAL_DIMENSION, T> > hessian_vals;
  this->evaluate_hessian_all(ips, hessian_vals);

  for ( j = 0; j < test_num; j++ )
  {
    T ddphi_sum = 0.0 ;
    for (i=0; i < nmodes; ++i)
    {
      cout << "  " << setw(12) << hessian_vals[j*nmodes+i].x(0) ;
      ddphi_sum += hessian_vals[j*nmodes+i].x(0) ;
    }
    cout << "  " << setw(12) << ddphi_sum << "\n";
  }

  hessian_vals.clear();

  return 0;
}


template <unsigned int ORDER, typename T>
inline
void JacobiSegmentFiniteElement<ORDER, T>::compute_vertices(std::vector<Point<SPATIAL_DIMENSION, T> >& vertices )
{
  std::vector<T> points(ORDER+1);
  weights.resize(ORDER+1);
  vertices.resize(ORDER+1);

  // eGaussLobattoLegendre
  polylib::zwglj(&points[0],&weights[0],ORDER+1,0.0,0.0);

  // eGaussGaussLegendre
  // polylib::zwgl(&points[0],&weights[0],ORDER+1,0.0,0.0);

  for (unsigned int i=0; i < ORDER+1; ++i)
    vertices[i][0] = points[i];

  points.clear();
}


template <unsigned int ORDER, typename T>
inline
void JacobiSegmentFiniteElement<ORDER, T>::reference_element_coords(std::vector<Point<SPATIAL_DIMENSION, T> >& nodes)
{
  const std::vector<Point<SPATIAL_DIMENSION, T> >& vertices = JacobiSegmentFiniteElement<ORDER, T>::get_vertices() ;

  nodes.resize(vertices.size());

  for (unsigned int i=0; i < vertices.size(); ++i)
    nodes[i][0] = vertices[i][0] ;
}



// Some useful typedefs
typedef JacobiSegmentFiniteElement<1, double> LinearJacobi1D ;
typedef JacobiSegmentFiniteElement<2, double> QuadraticJacobi1D ;
typedef JacobiSegmentFiniteElement<3, double> CubicJacobi1D ;
typedef JacobiSegmentFiniteElement<4, double> QuarticJacobi1D ;
typedef JacobiSegmentFiniteElement<5, double> QuinticJacobi1D ;



#endif // #ifndef __jacobi_fe_hh__
